#### Perform simulations using age-based yields

## Setup environment

setwd("~/research/coffee-dataverse/src")
source("simulate-yields.R")

do.origspec <- T
source("load.R", encoding="iso-8859-1")

output.suffix <- "-eq24-orig"

library(rstan)

## Load

param <- get.fit.params(output.suffix)
load(paste0("../data/combined-betagamma-arabica", output.suffix, ".RData"))

## Determine which municipalities grow only Arabica
muni.arabica <- unique(df.variety$Munic�pio[df.variety$arabica.portion == 1])
ii.arabica <- which(df$Munic�pio %in% muni.arabica)

sumstats <- read.csv("../data/sumstats-eq24-orig.csv")
names(sumstats)[1] <- "fullparam"
sumstats$param <- gsub("\\[\\d+\\]", "", gsub("_coeff\\[(\\d+)\\]", "_coeff\\1", sumstats$fullparam))

## Extract probability of die-off
deathprob <- mean(sumstats$mean[sumstats$region %in% muni.arabica & sumstats$param == 'death_part'])

## Extract other parameters
price <- mean(prices$price.arabica[prices$year >= 1980 & prices$year <= 2018])
cost <- mean(la$gamma[, 1, which(param == 'cost_harvest')])
delta <- mean(la$gamma[, 1, which(param == 'yield_uncertainty')])

## Simple models
model.annual <- function(ybar) {
    product <- rep(ybar, years)
    harvest <- rep(1, years)
    list(product=product, harvest=harvest)
}

model.byage <- function() {
    product <- yage.int
    harvest <- rep(1, years)
    list(product=product, harvest=harvest)
}

## NPV calculation
calc.npv <- function(result, discount) {
    sum((result$product * price - 1000 * result$harvest * cost) / ((1 + discount)^(0:(years-1))))
}

## Managed models

model.removal <- function(removeage) {
    result <- model.byage()
    product <- rep(result$product[1:removeage], floor(years / removeage) + 1)[1:years]
    list(product=product, harvest=result$harvest)
}

model.shocks <- function(result) {
    myrange <- result$product * delta
    product <- runif(years, result$product - myrange, result$product + myrange)
    list(product=product, harvest=result$harvest)
}

## Plotting model simulation
plot.result <- function(result, compare) {
    npv <- calc.npv(result, .05)

    if (length(compare) == 0)
        ggobj <- ggplot(data.frame(product=result$product, harvest=result$harvest), aes(x=1:100)) +
        geom_line(aes(y=product)) + geom_line(aes(y=harvest * 1000), linetype=2) +
        theme_bw()
    else
        ggobj <- ggplot(data.frame(product=result$product, compare, harvest=result$harvest), aes(x=1:100)) +
        geom_line(aes(y=product)) + geom_line(aes(y=harvest * 1000), linetype=2) +
        geom_line(aes(y=compare), linetype=3) +
        theme_bw()

    list(product=result$product, harvest=result$harvest, ggobj=ggobj, npv=npv)
}

## Highest NPV with removal at 33 years (20 if not interpolated)
## sapply(1:100, function(rm) calc.npv(model.removal(rm), .05))
## calc.npv(model.removal(33), .05)

## Plot simple models
plot.result(model.removal(33), model.byage()$product)
plot.result(model.shocks(model.removal(33)), model.removal(33)$product)

## Plot model requiring a Monte Carlo
plot.monte <- function(model, compare) {
    products <- matrix(NA, mciters, 100)
    harvests <- matrix(NA, mciters, 100)
    npvs <- c()
    for (ii in 1:mciters) {
        result <- model()
        products[ii,] <- result$product
        harvests[ii,] <- result$harvest
        npvs <- c(npvs, calc.npv(result, .05))
    }
    means <- colMeans(products)
    cilos <- as.numeric(sapply(1:100, function(ii) quantile(products[, ii], .25)))
    cihis <- as.numeric(sapply(1:100, function(ii) quantile(products[, ii], .75)))

    ggobj <- ggplot(data.frame(means, cilos, cihis, harvest=colMeans(harvests), compare),
           aes(x=1:100, y=means, ymin=cilos, ymax=cihis)) +
        geom_ribbon(alpha=.5) + geom_line(aes(y=means)) + geom_line(aes(y=harvest * 1000), linetype=2) +
        geom_line(aes(y=compare), linetype=3) +
        theme_bw()

    list(product=means, harvest=colMeans(harvests), ggobj=ggobj, npv=mean(npvs))
}

## Model die-off
model.deaths <- function(model) {
    product <- c()
    harvest <- c()
    while (length(product) < years) {
        alive <- rgeom(1, deathprob) + 1
        result <- model()
        product <- c(product, result$product[1:alive])
        harvest <- c(harvest, result$harvest[1:alive])
    }

    product <- product[1:years]
    harvest <- harvest[1:years]
    list(product=product, harvest=harvest)
}

compare <- plot.monte(function() model.shocks(model.removal(33)), model.removal(33)$product)$product
plot.monte(function() model.deaths(function() model.shocks(model.removal(33))), compare)

## Model harvest selection
model.harvestchoice <- function(result, minyield) {
    ## Harvest the entire plant or none of it
    product <- result$product
    product[result$product < minyield] <- 0
    harvest <- result$harvest
    harvest[result$product < minyield] <- 0

    list(product=product, harvest=harvest)
}

compare <- plot.monte(function() model.deaths(function() model.shocks(model.removal(33))), compare)$product
plot.monte(function() model.deaths(function() model.harvestchoice(model.shocks(model.removal(33)), cost/price)), compare)

## Model a random planting constraint
model.randomconst <- function(model, removeage) {
    product <- c()
    harvest <- c()
    while (length(product) < years) {
        result <- model()
        alive <- rgeom(1, deathprob) + 1
        if (alive > removeage) {
            probwait <- rep(.5, years)
            replant <- c(runif(length(probwait)) > probwait, T)
            replantage <- removeage + which(replant)[1] - 1
            if (alive >= replantage) {
                product <- c(product, result$product[1:replantage])
                harvest <- c(harvest, result$harvest[1:replantage])
            } else {
                product <- c(product, result$product[1:alive], rep(0, replantage - alive))
                harvest <- c(harvest, result$harvest[1:alive], rep(0, replantage - alive))
            }
        } else {
            probwait <- c(.5, years)
            replant <- c(runif(length(probwait)) > probwait, T)
            replantage <- alive + which(replant)[1] - 1

            product <- c(product, result$product[1:alive], rep(0, replantage - alive))
            harvest <- c(harvest, result$harvest[1:alive], rep(0, replantage - alive))
        }
    }

    product <- product[1:years]
    harvest <- harvest[1:years]
    list(product=product, harvest=harvest)
}

## Model a credit-based planting constraint
model.creditconst <- function(model, removeage) {
    product <- c()
    harvest <- c()
    while (length(product) < years) {
        result <- model()
        alive <- rgeom(1, deathprob) + 1
        if (alive > removeage) {
            probwait <- c(.8 - result$product[removeage:alive] / 1000, rep(.8, years))
            replant <- c(runif(length(probwait)) > probwait, T)
            replantage <- removeage + which(replant)[1] - 1
            if (alive >= replantage) {
                product <- c(product, result$product[1:replantage])
                harvest <- c(harvest, result$harvest[1:replantage])
            } else {
                product <- c(product, result$product[1:alive], rep(0, replantage - alive))
                harvest <- c(harvest, result$harvest[1:alive], rep(0, replantage - alive))
            }
        } else {
            probwait <- rep(.8, years)
            replant <- c(runif(length(probwait)) > probwait, T)
            replantage <- alive + which(replant)[1] - 1

            product <- c(product, result$product[1:alive], rep(0, replantage - alive))
            harvest <- c(harvest, result$harvest[1:alive], rep(0, replantage - alive))
        }
    }

    product <- product[1:years]
    harvest <- harvest[1:years]
    list(product=product, harvest=harvest)
}

## Standard graph saving logic
std.save <- function(filename, plotres) {
    ggobj2 <- plotres$ggobj + coord_cartesian(ylim=c(0, 1000)) + scale_x_continuous(expand=c(.01, 0)) +
        theme(axis.title.x=element_blank()) + ylab("Product, Harvest x 1000")
    ggsave(filename, ggobj2, width=4, height=2)

    print(plotres$npv)

    plotres$product
}

## Generate all graphs
setwd("../figures/sims")
compare <- std.save("model1-annual.pdf", plot.result(model.annual(ybar), c()))
compare <- std.save("model2-byage.pdf", plot.result(model.byage(), compare))
compare <- std.save("model3-removal.pdf", plot.result(model.removal(33), compare))
compare <- std.save("model4-shocks.pdf", plot.monte(function() model.shocks(model.removal(33)), compare))
compare <- std.save("model5-deaths.pdf", plot.monte(function() model.deaths(function() model.shocks(model.removal(33))), compare))
compare <- std.save("model6-harvestchoice.pdf", plot.monte(function() model.deaths(function() model.harvestchoice(model.shocks(model.removal(33)), cost/price)), compare))
compare <- std.save("model7-randomconst.pdf", plot.monte(function() model.randomconst(function() model.harvestchoice(model.shocks(model.byage()), cost/price), 33), compare))
std.save("model8-creditconst.pdf", plot.monte(function() model.creditconst(function() model.harvestchoice(model.shocks(model.byage()), cost/price), 33), compare))

## Consider a range of death rates
## For each, extract (1) ratio of mean in 32 and 33 under deaths
## NPV after deaths and after harvestchoice

mciters <- 10000
results <- data.frame(deathprob=c(), npv.death=c(), npv.harvest=c(), level.ratio=c(), product.ratio=c())
for (deathprob in exp(seq(-6, 0, length.out=50))) {
    after.death <- plot.monte(function() model.deaths(function() model.shocks(model.removal(33))), rep(0, 100))
    after.harvest <- plot.monte(function() model.deaths(function() model.harvestchoice(model.shocks(model.removal(33)), cost/price)), rep(0, 100))

    results <- rbind(results, data.frame(deathprob, npv.death=after.death$npv, npv.harvest=after.harvest$npv, level.ratio=after.death$product[34] / after.death$product[33], harvest=mean(after.harvest$harvest)))
    print(results[nrow(results),])
}
results$level.ratio[is.nan(results$level.ratio) | results$level.ratio > 1] <- 1

labels <- c("NPV with full harvest", "NPV after harvest selection", "Removal ratio", "Harvest ratio")

## Produce graph

disps <- data.frame(deathprob=rep(results$deathprob, 4), values=c(results$npv.death, results$npv.harvest, lowess(results$level.ratio, f=.2)$y, results$harvest), panel=c(rep("Net Present Value", 2*nrow(results)), rep("Characteristic Ratio", 2*nrow(results))), group=rep(labels, each=nrow(results)))

disps2 <- rbind(subset(disps, group != 'Removal ratio'),
                data.frame(deathprob=exp(seq(-6, 0, length.out=50)), values=1, panel="Harvest Ratio", group="Full harvest"))
disps2$group <- as.character(disps2$group)
disps2$group[disps2$group == "Harvest ratio"] <- "Harvest selection ratio"
disps2$group <- factor(disps2$group, levels=c("NPV with full harvest", "NPV after harvest selection", "Full harvest", "Harvest selection ratio"))
disps2$panel <- as.character(disps2$panel)
disps2$panel[disps2$panel == "Characteristic Ratio"] <- "Harvest Ratio"
disps2$panel <- factor(disps2$panel, levels=c("Net Present Value", "Harvest Ratio"))

ggplot(disps2, aes(x=deathprob, y=values, colour=group, linetype=group)) +
    facet_grid(panel ~ ., scales="free") +
    scale_colour_discrete(name="Metric") +
    scale_linetype_manual(name="Metric", values=c('solid', 'dashed', 'solid', 'dashed')) +
    xlab("Probability of die-off") + scale_x_log10() +
    geom_line() + theme_bw() + theme(axis.title.y=element_blank())
ggsave("dierange.pdf", width=6.5, height=4.5)
